rm(list=ls())
library(epicalc)
zap()
#Spatial spec correct, omitted variable in X\betas (correlated with \rho)

sar.f <- function(es, x, y, w){
    rho <- es[1]
    s2 <- es[2]^2
    if (s2<0) return(-1e20)
    k <- ncol(x)
    b <- as.matrix(es[3:(k+2)])
    llik <- 0
    nd <- nrow(w)
    n <- nrow(x)
    d <- n/nd
    y <- as.matrix(y)
    
    for(i in 1:d){
        A <- diag(1, nd, nd)-rho*w[,,i]
        dA <- det(A)
        if (dA <=0){
            llik <-  -1e20 
            break
        } else {
            aux <- log(dA) -nd/2 * log(s2)- (1/(2*s2)) * t(A %*% y[((i-1)*nd+1):((i-1)*nd+nd),] - as.matrix(x[((i-1)*nd+1):((i-1)*nd+nd),])%*%b) %*% (A %*% y[((i-1)*nd+1):((i-1)*nd+nd),]- as.matrix(x[((i-1)*nd+1):((i-1)*nd+nd),])%*%b)
            llik <- llik + aux
        }
    }
    return(llik)
}

sar <- function(x, y, w, stval,vars=c(""), routit=100000){
    mod <- optim(stval, sar.f, hessian = TRUE, method = "BFGS", 
            control = list(fnscale = -1, REPORT = 10, 
            reltol = 1e-16, trace = 1, 
            maxit = routit), x=x, y=y, w=w)
    vacov <- as.matrix(solve(-1*mod$hessian, tol=1e-24))
    se <- as.matrix(sqrt(diag(vacov)))
    zstat <- as.matrix(mod$par / se)
    pvalues <- round(2*pnorm(abs(zstat),lower.tail=FALSE),4)
    
    cat("results in order: rho, s, b \n")
    cat("parameter estimates: ", mod$par, "\n")
    cat("standard error: ", se, "\n")
    cat("z-statistics: ", zstat, "\n")
    cat("p-values: ", pvalues, "\n")
    list(llik=mod$value, par=mod$par, se=se, zstat=zstat, 
        pvalues = pvalues, 
        hessian=mod$hessian, vacov=vacov
        )
}



genY <- function(nd, d, ld, sig2, rho, b0, b1, b2){ #Generate data from correct DGP

    wall.ld <- array(NA, c(nd, nd, d))
    wraw <- matrix(c(rep(.5/(nd-2),nd)),nd, nd, byrow=TRUE)
    ldrow <- rep(1/(nd-1), nd)
    for (i in 1:d){
        wsub <- wraw
        wsub[ld[i],] <- ldrow
        wsub[,ld[i]] <- rep(.5, nd)
        diag(wsub) <- 0
        wall.ld[,,i] <- wsub
    }    
        
    ytrue <- rep(NA, d*nd)
    xldtrue <- rep(NA, d*nd)
    xother <- rep(NA, d*nd)
    for (j in 1:d){
        xld <- rep(0,nd)
        xld[ld[j]] <- 1
        iA <- solve(diag(1,nd)-rho * wall.ld[,,j])
        xotheraux <- rnorm(nd)
        yaux <- iA %*% (b0 + b1 * xld+ b2 * xotheraux) + iA%*%rnorm(nd, sd=sqrt(sig2))
        xldtrue[(1+(j-1)*nd):(nd+(j-1)*nd)] <- xld
        ytrue[(1+(j-1)*nd):(nd+(j-1)*nd)]<- yaux
        xother[(1+(j-1)*nd):(nd+(j-1)*nd)]<- xotheraux
    }
    return(list(ytrue=ytrue, xldtrue=xldtrue, xother=xother, wall.ld=wall.ld))
}


#repetitions
r <- 1000
    
#specify parameters
sig2 <- 1.0
rho <- -.5
b0 <- .5
b1 <- 1
b2 <- 2

nd<-10
d<-100
set.seed(021175)

mod.true <- vector("list", r)
mod.lack <- vector("list", r)
i <- 1
while (i <= r){
    #vector of lead donors
    dn <- 1:nd
    ld <- sample(dn, d, TRUE)
    
    mcdat <- genY(nd, d, ld, sig2, rho, b0, b1, b2)
    xfull <- cbind(1, mcdat$xldtrue, mcdat$xother)
    xlack <- cbind(1, mcdat$xother)
    y <- mcdat$ytrue
    stvfull <- runif((ncol(xfull) + 2))*.1
    stvlack <- runif((ncol(xlack) + 2))*.1
    mod.true[[i]] <- try(sar(xfull, y, mcdat$wall.ld, stval=stvfull, vars=c("const", "ld", "other")), TRUE)
    ect <- 1 #initiate error counter
    while(class(mod.true[[i]]) == "try-error"){
            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
            stvfull <- runif((ncol(xfull) + 2))*.1
            mod.true[[i]] <- try(sar(xfull, y, mcdat$wall.ld, stval=stvfull, vars=c("const", "ld", "other")), TRUE)
            ect <- ect+1
            if (ect == 10) {
                i <-  r+1
                cat("\n\n\n\n ROUTINE FAILED BC OF TRUE MODEL\n")
                break
            }
    }    
    mod.lack[[i]] <- try(sar(xlack, y, mcdat$wall.ld, stval=stvlack, vars=c("const", "other")), TRUE)
    ect <- 1 #initiate error counter
    while(class(mod.lack[[i]]) == "try-error"){
            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
            stvlack <- runif((ncol(xlack) + 2))*.1
            mod.lack[[i]] <- try(sar(xlack, y, mcdat$wall.ld, stval=stvlack, vars=c("const", "other")), TRUE)
            ect <- ect+1
            if (ect == 10) {
                i <-  r+1
                cat("\n\n\n\n ROUTINE FAILED BC OF FALSE MODEL\n")
                break
            }
    }
    cat("\n\n finished repetition ", i, " \n\n")
    i <- i + 1
}

getpars <- function(i, model){
    model[[i]]$par
}
getllik <- function(i, model){
    model[[i]]$llik
}
pars.true <- unlist(sapply(1:1000, getpars, mod.true), recursive=FALSE)
pars.true <- matrix(pars.true, r, length(stvfull), byrow=TRUE)

pars.lack <- unlist(sapply(1:1000, getpars, mod.lack), recursive=FALSE)
pars.lack <- matrix(pars.lack, r, length(stvlack), byrow=TRUE)

llik.true <- unlist(sapply(1:r, getllik, mod.true), recursive=FALSE)
llik.lack <- unlist(sapply(1:r, getllik, mod.lack), recursive=FALSE)

quantile(pars.true[,1], c(.025,.5,.975))
quantile(pars.lack[,1], c(.025,.5,.975))

bias.true <- mean(pars.true[,1]-rho)
bias.lack <- mean(pars.lack[,1]-rho)

mean(llik.true)
mean(llik.lack)


#true model has better fit:
#ttest of likelihood
t.test(llik.true,llik.lack,paired=TRUE)

#llratio test of mean llik
dgf=1
r <- -2 * (mean(llik.lack)-mean(llik.true))
p <- 1-pchisq(r, dgf)
cat("degrees of freedom: ", dgf, "\n")
cat("likelihood ratio statistic: ", r,"\n")
cat("p value: ", p, " (under H0 that the restricted model is correct)\n\n")
# number of lliks greater for full model
check<- llik.true>llik.lack
sum(check)

setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/R&R.IO")
#save(mod.lack, file="mcld.omit.out")
#save(mod.true, file="mcld.full.out")

##################
## Show that variance does not matter for retrieved values of \rho
##################

#W matrix -- even
#array

genY.var <- function(nd, d, ld, sig2, rho, b0, b2, lamb){ #Generate data from correct DGP

    wall.ev <- array(NA, c(nd, nd, d))
    wraw <- matrix(c(rep(1/(nd-1),nd)),nd, nd, byrow=TRUE)
    diag(wraw) <- 0
    wall.ev[,,1:d] <- wraw        
    ytrue <- rep(NA, d*nd)
    xother <- rep(NA, d*nd)
    befftrue <- rep(NA, d*nd)
    for (j in 1:d){
        beff <- b0+rexp(1,lamb)
        iA <- solve(diag(1,nd)-rho * wall.ev[,,j])
        xotheraux <- rnorm(nd)
        yaux <- iA %*% (beff +  b2 * xotheraux) + iA%*%rnorm(nd, sd=sqrt(sig2))
        ytrue[(1+(j-1)*nd):(nd+(j-1)*nd)]<- yaux
        befftrue[(1+(j-1)*nd):(nd+(j-1)*nd)]<- beff
        xother[(1+(j-1)*nd):(nd+(j-1)*nd)]<- xotheraux
    }
    return(list(ytrue=ytrue, befftrue=befftrue, xother=xother, wall.ev=wall.ev))
}


#repetitions
r <- 1000
    
#specify parameters
sig2 <- 1.0
rho <- .5
b0 <- .5
b2 <- 2
lamb.l <- 2
lamb.h <- .5

nd<-10
d<-100
set.seed(021175)

mod.l <- vector("list", r)
mod.h <- vector("list", r)
i <- 1

while (i <= r){
    #vector of lead donors
    dn <- 1:nd
    ld <- sample(dn, d, TRUE)
    
    mcdat.l <- genY.var(nd, d, ld, sig2, rho, b0, b2, lamb.l)
    mcdat.h <- genY.var(nd, d, ld, sig2, rho, b0, b2, lamb.h)
    Xsim.l <- cbind(1, mcdat.l$xother)
    Xsim.h <- cbind(1, mcdat.h$xother)
    y.l <- mcdat.l$ytrue
    y.h <- mcdat.h$ytrue
    stv <- runif((ncol(Xsim.l) + 2))*.1
    mod.l[[i]] <- try(sar(Xsim.l, y.l, mcdat.l$wall.ev, stval=stv, vars=c("const", "other")), TRUE)
    ect <- 1 #initiate error counter
    while(class(mod.l[[i]]) == "try-error"){
            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
            stv <- runif((ncol(Xsim.l) + 2))*.1
            mod.l[[i]] <- try(sar(Xsim.l, y.l, mcdat.l$wall.ev, stval=stv, vars=c("const", "other")), TRUE)
            ect <- ect+1
            if (ect == 10) {
                i <-  r+1
                cat("\n\n\n\n ROUTINE FAILED BC OF LOW VAR MODEL\n")
                break
            }
    }    
    mod.h[[i]] <- try(sar(Xsim.h, y.h, mcdat.h$wall.ev, stval=stv, vars=c("const", "other")), TRUE)
    ect <- 1 #initiate error counter
    while(class(mod.h[[i]]) == "try-error"){
            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
            stv <- runif((ncol(Xsim.h) + 2))*.1
            mod.h[[i]] <- try(sar(Xsim.h, y.h, mcdat.h$wall.ev, stval=stv, vars=c("const", "other")), TRUE)
            ect <- ect+1
            if (ect == 10) {
                i <-  r+1
                cat("\n\n\n\n ROUTINE FAILED BC OF HIGH VAR MODEL\n")
                break
            }
    }
    cat("\n\n finished repetition ", i, " \n\n")
    i <- i + 1
}

setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/R&R.IO")
save(mod.l, file="mcvar.low.out")
save(mod.h, file="mcvar.high.out")


getpars <- function(i, model){
    model[[i]]$par
}
pars.l <- unlist(sapply(1:1000, getpars, mod.l), recursive=FALSE)
pars.l <- matrix(pars.l, r, length(stv), byrow=TRUE)

pars.h <- unlist(sapply(1:1000, getpars, mod.h), recursive=FALSE)
pars.h <- matrix(pars.h, r, length(stv), byrow=TRUE)





############################################################
### Demonstrate simple scaling effect affects \rho estimates
genY.scale <- function(nd, d, ld, sig2, rho, b0, b2){ #Generate data from correct DGP

    wall.ev <- array(NA, c(nd, nd, d))
    wraw <- matrix(c(rep(1/(nd-1),nd)),nd, nd, byrow=TRUE)
    diag(wraw) <- 0
    wall.ev[,,1:d] <- wraw        
    ylow <- rep(NA, d*nd)
    yhigh <- rep(NA, d*nd)
    xother <- rep(NA, d*nd)
    for (j in 1:d){
        iA <- solve(diag(1,nd)-rho * wall.ev[,,j])
        xotheraux <- rnorm(nd)
        yaux <- iA %*% (b0 +  b2 * xotheraux) + iA%*%rnorm(nd, sd=sqrt(sig2))
        ylow[(1+(j-1)*nd):(nd+(j-1)*nd)]<- yaux
        yhigh[(1+(j-1)*nd):(nd+(j-1)*nd)]<- yaux*1000
        xother[(1+(j-1)*nd):(nd+(j-1)*nd)]<- xotheraux
    }
    return(list(yhigh=yhigh, ylow=ylow, xother=xother, wall.ev=wall.ev))
}


#repetitions
r <- 1000
    
#specify parameters
sig2 <- 1.0
rho <- .5
b0 <- .5
b2 <- 2

nd<-10
d<-100
set.seed(021175)

mod.l <- vector("list", r)
mod.h <- vector("list", r)
i <- 1

while (i <= r){
    #vector of lead donors
    dn <- 1:nd
    ld <- sample(dn, d, TRUE)
    
    mcdat <- genY.scale(nd, d, ld, sig2, rho, b0, b2)
    Xsim <- cbind(1, mcdat$xother)
    y.l <- mcdat$ylow
    y.h <- mcdat$yhigh
    stv <- runif((ncol(Xsim) + 2))*.1
    mod.l[[i]] <- try(sar(Xsim, y.l, mcdat$wall.ev, stval=stv, vars=c("const", "other")), TRUE)
    ect <- 1 #initiate error counter
    while(class(mod.l[[i]]) == "try-error"){
            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
            stv <- runif((ncol(Xsim.l) + 2))*.1
            mod.l[[i]] <- try(sar(Xsim, y.l, mcdat$wall.ev, stval=stv, vars=c("const", "other")), TRUE)
            ect <- ect+1
            if (ect == 10) {
                i <-  r+1
                cat("\n\n\n\n ROUTINE FAILED BC OF LOW SCALE MODEL\n")
                break
            }
    }    
    mod.h[[i]] <- try(sar(Xsim, y.h, mcdat$wall.ev, stval=stv, vars=c("const", "other")), TRUE)
    ect <- 1 #initiate error counter
    while(class(mod.h[[i]]) == "try-error"){
            cat("\n\n restarting repetition ",i,", ",ect,". attempt \n\n")
            stv <- runif((ncol(Xsim) + 2))*.1
            mod.h[[i]] <- try(sar(Xsim, y.h, mcdat$wall.ev, stval=stv, vars=c("const", "other")), TRUE)
            ect <- ect+1
            if (ect == 10) {
                i <-  r+1
                cat("\n\n\n\n ROUTINE FAILED BC OF HIGH VAR MODEL\n")
                break
            }
    }
    cat("\n\n finished repetition ", i, " \n\n")
    i <- i + 1
}

setwd("C:/Users/Martin Steinwand/Documents/Project Lead Donor/Analysis/R&R.IO")
#save(mod.l, file="mcscale.low.out")
#save(mod.h, file="mcscale.high.out")


getpars <- function(i, model){
    model[[i]]$par
}
pars.l <- unlist(sapply(1:1000, getpars, mod.l), recursive=FALSE)
pars.l <- matrix(pars.l, r, length(stv), byrow=TRUE)

pars.h <- unlist(sapply(1:1000, getpars, mod.h), recursive=FALSE)
pars.h <- matrix(pars.h, r, length(stv), byrow=TRUE)
